{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<script>\n",
       "code_show=true; \n",
       "function code_toggle() {\n",
       " if (code_show){\n",
       " $('div.input').hide();\n",
       " } else {\n",
       " $('div.input').show();\n",
       " }\n",
       " code_show = !code_show\n",
       "} \n",
       "$( document ).ready(code_toggle);\n",
       "</script>\n",
       "<form action=\"javascript:code_toggle()\"><input type=\"submit\" value=\"Click here to toggle on/off the raw code for QISKit exercises.\"></form>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import HTML\n",
    "\n",
    "HTML('''<script>\n",
    "code_show=true; \n",
    "function code_toggle() {\n",
    " if (code_show){\n",
    " $('div.input').hide();\n",
    " } else {\n",
    " $('div.input').show();\n",
    " }\n",
    " code_show = !code_show\n",
    "} \n",
    "$( document ).ready(code_toggle);\n",
    "</script>\n",
    "<form action=\"javascript:code_toggle()\"><input type=\"submit\" value=\"Click here to toggle on/off the raw code for QISKit exercises.\"></form>''')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 8. Shor's algorithm\n",
    "\n",
    "\n",
    "Shor's factoring algorithm [1] is the most well-known example of a quantum algorithm outperforming the best known classical algorithms. The algorithm allows to factor a number $N_0$ which is the product of two prime numbers $N_0 = p \\cdot q$, in polynomial time. This is possible thanks to a theorem in number theory which turns the problem of finding factors into the problem of finding the period of a periodic function. Using the quantum Fourier transform (QFT), one can then find the period of this function with high probability and complete the factorization.\n",
    "\n",
    "The problem of factoring numbers has been studied for centuries and no efficient (polynomial time) algorithm has ever been found. The difficulty of factoring numbers is the basis on which the most widespread encryption standard, the RSA, is founded. Therefore, an algorithm which is capable of factoring numbers efficiently could have a huge impact on the security of electronic based interactions. The discovery of this algorithm by Peter Shor in 1994, led to an explosion of the field of quantum computation because of its important application. \n",
    "\n",
    "However, experimental implementation of the algorithm still remains a challenge because of the errors introduced by the large number of physical qubits and gates required to execute the algorithm. Proof of principle demonstration of Shor's factoring algorithm to factor the smallest number $N_0=15$ using five superconducting qubits have been done for setups like NMR [2], trapped ions [3], photons [4-6], photonic chips [7] and superconducting qubits [8,9]. These experiments show how complex is to implement the actual algorithm, as even for the smallest number which allows the algorithm to be run ($N_0=15$), extreme simplifications must be made to run the algorithm [10]. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.1 Mathematical preliminaries\n",
    "\n",
    "\n",
    "\n",
    "### 8.1.1 Modular arithmetic\n",
    "When division between integers $a$ (the dividend) and $b$ (the divisor) is taught for the first time, one is usually introduced to the idea of quotient $q$ and remainder $r$. The quotient is the number of times the dividend $a$ contains the divisor $b$, the remainder is any left-over that added to the product of the quotient by the divisor returns the dividend, $r + q\\times b = a$ or $r = a - q\\times b$.\n",
    "\n",
    "#### Example:\n",
    "Consider the case $a=40$ and $b=30$, let us find the quotient and the remainder of the division of $a$ by $b$. $40/30 \\rightarrow q = 1, \\, r = 13$,\n",
    "\t\t \n",
    "Therefore, the quotient is $q=1$ and the remainder is $r = 13$. \n",
    "\t\n",
    "\n",
    "\n",
    "By looking at the division between integers in this way, one can easily understand the fundamentals of modular arithmetic. Modular arithmetic is a set of rules for handling operations between integer numbers. In modular arithmetic, integers multiple of a fixed number $N_0$ are considered equivalent. Thus, by selecting an integer $N_0$, the set of all integers is restricted to the integers in the interval $\\left[ 0 , N_0-1 \\right]$. \n",
    "To visualize this, one can think about a twelve-hours clock. In a clock, the only integers allowed lie in the interval $\\left[ 0 , 11 \\right]$, once this interval is exceeded one goes back at the beginning of the interval. Thus, $9+5 = 2 \\, \\left( mod \\, 12 \\right)$.\n",
    "More formally, two integers $a$ and $b$ are congruent modulo $N_0$ if their difference $a-b$ is an integer multiple of $N_0$. That is, $a-b = r \\cdot N_0$, where $r$ is an integer, then $a \\equiv b  \\, \\left( mod \\, N_0 \\right)$. \n",
    "Equivalently, if $a \\equiv b  \\, \\left( mod \\, N_0 \\right)$, then $\\frac{a-b}{N_0}$ has zero remainder. Or, $b  \\, \\left( mod \\, N_0 \\right)$ is the remainder of the division $\\frac{b}{N_0}=a$. Thus, both numbers have the same remainder when divided by $N_0$.\n",
    "\n",
    "#### Example:\n",
    "Equivalence of two integers modulo $N_0$:\n",
    "$38 \\equiv 14  \\, \\left( mod \\, 12 \\right)$\n",
    "$38 - 14  = 24 = 2 \\cdot 12 $\n",
    "\n",
    "Both $\\frac{38}{12}$ and $\\frac{14}{12}$ have remainder $2$.\n",
    "\n",
    "\n",
    "\n",
    "### 8.1.2 Continued Fraction Algorithm\n",
    "\n",
    "The continued fraction algorithm is used to reduce a fraction $\\frac{m}{n}$ to another fraction $\\frac{u}{t}$. In general, the algorithm allows to rewrite any irrational number as a finite/infinite sum of an integer part plus the reciprocal of a number.\n",
    "So, consider the fraction $\\frac{m}{n}$, its continued fraction expansion is\n",
    "\n",
    "$$ \\frac{m}{n} = a_0 + \\frac{1}{a_1 + \\frac{1}{a_2+ \\frac{1}{a_3+\\frac{1}{1+\\frac{1}{a_4+\\frac{1}{a_5+\\frac{1}{...}}}}}}} $$\n",
    "\n",
    "Using the integers $a_i$, one can rewrite the fraction $\\frac{m}{n}$ as $\\frac{u}{t}$. To find $u_i$ and $t_i$ one can use the following formulas:\n",
    "\n",
    "$u_0 = a_0$, $u_1=1+a_0a_1$, ..., $u_n = a_n d_n + d_{n-2}$\n",
    "$t_0=1$, $t_1 = a_1$, ..., $t_n =a_n t_{n−1} +t_{n−2}$ \n",
    "\n",
    "Which then give the possible approximations for $\\frac{m}{n}$:\n",
    "\n",
    "$\\frac{m}{n} \\approx \\frac{u_0}{t_0}$, $\\frac{m}{n} \\approx \\frac{u_1}{t_1}$, ..., $\\frac{m}{n} \\approx \\frac{u_n}{t_n}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.2 Factoring and order finding \n",
    " \n",
    " \n",
    "The possibility of running an efficient algorithm for factoring a product of two prime number arises from: i. the connection between factoring and order finding, ii. the ability of quantum computers to deal efficiently with periodic functions.\n",
    "In this Section, the connection between the problem of factoring and the problem of finding the period of a periodic function, also called $\\textit{order finding}$ is explained.\n",
    "It is important to note that this equivalency holds in the realm of modular arithmetic.\n",
    "Let us start from the concept of order of a number. Given an integer $a$, the order of $a$ is the smallest integer number $r$ for which the following condition holds\n",
    "\n",
    "\\begin{equation}\n",
    "a^r \\text{mod} N_0 = 1 \\tag{1}\n",
    "\\end{equation}\n",
    " \n",
    "\n",
    "\n",
    "#### Example:\n",
    "Consider the case $a=2$ and $N_0=21$, let us find the order of $a$ by trying different exponents $r$ sequentially until \twe find the one for which condition (1) is satisfied:\n",
    "\t$2^1 (\\text{mod}  21) = 2$, \n",
    "\t$2^2 (\\text{mod}  21) = 4$, \n",
    "\t$2^3 (\\text{mod} 21) = 8$, \n",
    "\t$2^4 (\\text{mod} 21) = 16$,\n",
    "\t$2^5 (\\text{mod} 21) = 11$, \n",
    "\t$2^6 (\\text{mod} 21) = 1$.\n",
    "\n",
    "Therefore, the order of $a=2$ is $r=6$. \n",
    "\t\n",
    "\n",
    "\n",
    "Now that we have clarified what is the order of a number, let us explore the connection between order finding and factoring. \n",
    "Start by rewriting Eq. (1) as $a^r - 1  \\text{mod} \\, N_0 = 0$, which means that $a^r - 1$ is a multiple of $N_0$. Assuming that $r$ is even, we can write $ a^r - 1  \\text{mod} \\, N_0  = \\left( a^{\\frac{r}{2}} - 1\\right)\\left( a^{\\frac{r}{2}} + 1\\right)  \\text{mod} \\, N_0 $. Finally, we check if either $a^{\\frac{r}{2}} - 1$ or $a^{\\frac{r}{2}} + 1$ shares any factors with $N_0$. This can be checked tacking the greatest common divisor (gcd) of them, which outputs the biggest number that divides both inputs. If the output of $\\text{gcd}\\left( a^{\\frac{r}{2}} \\pm 1,N_0 \\right)$ is greater than one, than that's one of the factors. The other one being $N_0/ \\text{gcd}\\left( a^{\\frac{r}{2}} \\pm 1,N_0 \\right)$.\n",
    "\n",
    "#### Example:\n",
    "Again consider the case $a=2$ and $N_0=21$, in the previous example we found that the order is $r = 6$. Let us now show how to find the two prime factors whose product gives $21$. \n",
    "$2^6 - 1 (\\text{mod} \\, 21) = \\left( 2^3 +1 \\right)\\left( 2^3 -1 \\right) (\\text{mod} \\, 21) $.\n",
    "We need to check if the two numbers $ \\left( 2^3 +1 \\right) = 9$ and $ \\left( 2^3 -1 \\right)=7$ have any factors in common with $21$.\n",
    "Let's start by checking $ \\left( 2^3 +1 \\right)$: $\\text{gcd}\\left( 9, 21\\right) = 3$.\n",
    "We obtained a $gcd \\neq 1$, which means that there is a common factor between $a^{\\frac{r}{2}} + 1$ and $N_0$. We can now calculate the first factor of $21$.\n",
    "$\\text{factor}_1 = \\frac{21}{\\text{gcd}\\left( 9, 21\\right)} = \\frac{21}{3} = 7$.\n",
    "The second factor follows immediately from $\\text{factor}_2 = \\frac{21}{\\text{factor}_1} = \\frac{21}{7} = 3$.\n",
    "\n",
    "Therefore, we find the two prime factors whose product is $N_0=21$, $\\text{factor}_1 = 7$ and $\\text{factor}_2 = 3$. Thus, the order of $a=2$ is $r=6$. \n",
    "\t\n",
    "\n",
    "\n",
    "It's important to notice that several assumptions are needed to translate the problem of factoring to the one of order finding. These are \"weak points\" of the algorithm, which fails to produce the correct result each time one of these assumptions is violated.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.3 Quantum Fourier transform\n",
    "\n",
    "The Quantum Fourier Transform (QFT) is the heart of Shor's factoring algorithm. Similarly to the discrete Fourier transform, which allows to compute the Fourier transform of a complex number, the QFT allows to compute the Fourier transform of a quantum state. This means that we can rewrite a generic state of a qubit register as a superposition of all the possible basis state vectors of the register, with a certain phase.\n",
    "\n",
    "The QFT of one of the basis states $\\lvert i \\rangle$ of an $n$ qubit register is\n",
    "\n",
    "\\begin{equation}\n",
    "QFT \\lvert j \\rangle = \\frac{1}{2^{n/2}} \\sum_{k=0}^{2^{n}-1} e^{\\frac{2 \\pi i}{2^n} j k} \\lvert k \\rangle\n",
    "\\end{equation}\n",
    "\n",
    "The QFT on a generic state can then be derived from this definition\n",
    "\n",
    "\\begin{equation}\n",
    "QFT \\sum_{j=0}^{2^n-1} x_j \\lvert j \\rangle = \\sum_{k=0}^{2^{n}-1}  \\sum_{l=0}^{2^{n}-1} x_l e^{\\frac{2 \\pi i}{2^n} l k} \\lvert k \\rangle\n",
    "\\end{equation}\n",
    "\n",
    "By writing the basis state $\\lvert j \\rangle$ in binary representation we adopt the following convention: $\\lvert j_1 j_2 ... j_n \\rangle$ correspond to $j= j_1 2^{n-1}+j_2 2^{n-2} ... j_n 2^{0}$. Following the same convention, we can also write numbers smaller than one as $0.j_1 j_2 ... j_n$  meaning $j_1/2^{0} + j_2/2^{1} +...+ j_n/2^{n}$ a  it is possible to rewrite the QFT in a very simple form\n",
    "\n",
    "\\begin{aligned}\n",
    "QFT \\lvert j \\rangle & = \\frac{1}{2^{n/2}} \\sum_{k=0}^{2^{n}-1} e^{2 \\pi i jk / 2^n} \\lvert k \\rangle \\\\\n",
    "& = \\frac{1}{2^{n/2}} \\sum_{k_1=0}^{1}...\\sum_{k_n=0}^{1} e^{2 \\pi i \\left(\\sum_{l=1}^{n} k_l 2^{-l}\\right) j} \\lvert k_1 ... k_n \\rangle \\\\\n",
    "& = \\frac{1}{2^{n/2}}\\sum_{k_1=0}^{1}...\\sum_{k_n=0}^{1} \\bigotimes_{l=0}^n e^{2 \\pi i j k_l2^{-l} } \\lvert k_l \\rangle \\\\\n",
    "& = \\frac{1}{2^{n/2}} \\bigotimes_{l=1}^n  \\left(\\lvert 0 \\rangle + e^{2 \\pi i j 2^{-l} } \\lvert 1 \\rangle \\right) \\\\\n",
    "& = \\frac{1}{2^{n/2}} \\left( \\lvert 0 \\rangle + e^{2 \\pi i[0.j_n]} \\lvert 1 \\rangle \\right) \\otimes \\left( \\lvert 0 \\rangle + e^{2 \\pi i[0.j_{n-1} j_{n}]} \\lvert 1 \\rangle \\right) \\otimes ...\\otimes  \\left(\\lvert 0 \\rangle + e^{2 \\pi i[0.j_1j_2...j_{n-1}j_n]} \\lvert 1 \\rangle \\right)\n",
    "\\end{aligned}\n",
    "\n",
    "By writing the QFT in this way, it is easy to find a circuit which implements it. The QFT turns out to be composed of very simple gates such as Hadamard and control-rotations around the z-axis between pairs of qubits $CR_k$.\n",
    "\n",
    "$$ CR_k =\n",
    "\\begin{pmatrix} \n",
    "1 & 0 & 0 & 0 \\\\\n",
    "0 & 1 & 0 & 0 \\\\ \n",
    "0 & 0 & 1 & 0 \\\\ \n",
    "0 & 0 & 0 & e^{2 \\pi i/ 2^k}\n",
    "\\end{pmatrix}$$\n",
    "\n",
    "<img src=\"figures/8/QFT1.jpeg\"  width=\"900\">\n",
    "$$\\text{1. Quantum circuit for the Quantum Fourier Transform.}$$\n",
    "\n",
    "\n",
    "#### Example: QFT on two qubits\n",
    "Let us calculate the QFT in the case of $n=2$ qubits\n",
    "\n",
    "$$QFT\\lvert j_1j_2\\rangle = \\frac{1}{\\sqrt{8}} \\left(\\lvert 0 \\rangle + e^{2 \\pi i[0.j_2]} \\lvert 1 \\rangle\\right) \\otimes \\left(\\lvert 0 \\rangle + e^{2 \\pi i[0.j_1 j_2]} \\lvert 1 \\rangle \\right) $$\n",
    "\n",
    "The steps to creating the circuit for $\\lvert k_1k_2\\rangle = QFT\\lvert j_1j_2\\rangle$ would be:\n",
    "\n",
    "<ol>\n",
    "<li> \n",
    "Apply an Hadamard to $\\lvert j_2 \\rangle$, giving the state $\\frac{1}{\\sqrt{2}}\\left(\\lvert 0 \\rangle + e^{2 \\pi i 0.j_2} \\lvert 1\\rangle\\right) = \\frac{1}{\\sqrt{2}}\\left(\\lvert 0 \\rangle + (-1)^{j_2} \\lvert 1 \\rangle\\right)$\n",
    "</li>\n",
    "\n",
    "<li> \n",
    "Apply a $CR_2$ from $\\lvert j_1 \\rangle$ to $\\lvert j_2 \\rangle$ and an Hadamard to $\\lvert j_1 \\rangle$, giving the state $\\frac{1}{\\sqrt{2}}\\left(\\lvert0\\rangle + e^{2 \\pi i[0.j_1j_2]} \\lvert1\\rangle\\right)$\n",
    "</li>\n",
    "\n",
    "<li> \n",
    "Measure the qubits in reverse order, that is $k_2$ corresponds to the first qubit, $k_1$ to the second qubit.\n",
    "</li>\n",
    "</ol>\n",
    "\n",
    "<img src=\"figures/8/QFT_ex1.jpeg\"  width=\"250\">\n",
    "$$\\text{2. Quantum Fourier Transform for two qubits.}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# QISKIT: QFT code \n",
    "\n",
    "# calculate qft of two qubit as the example above"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color:blue\"> QISKit: implement the quantum Fourier transform </span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <span style=\"color:blue\"> 1) Fourier transform of two qubits in the $\\lvert 00 \\rangle $ state </span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "List of gates:\n",
      "h\n",
      "cu1\n",
      "h\n",
      "measure\n",
      "measure\n",
      "\n",
      "The measured outcomes of the circuits are: {'00': 243, '01': 241, '10': 272, '11': 268}\n"
     ]
    }
   ],
   "source": [
    "from initialize import *\n",
    "\n",
    "def qft(circ, q, n):\n",
    "    \"\"\"n-qubit QFT on q in circ.\"\"\"\n",
    "    for j in range(n):\n",
    "        for k in range(j):\n",
    "            circ.cu1(3.14/float(2**(j-k)), q[j], q[k])\n",
    "        circ.h(q[j])\n",
    "\n",
    "#initialize quantum program\n",
    "my_alg = initialize(circuit_name = 'fourier', qubit_number=2, bit_number=2, backend = 'local_qasm_simulator', shots = 1024)\n",
    "\n",
    "#add gates to the circuit\n",
    "qft(my_alg.q_circuit, my_alg.q_reg, 2) # qft of input state\n",
    "my_alg.q_circuit.measure(my_alg.q_reg[0], my_alg.c_reg[0]) # measures first qubit\n",
    "my_alg.q_circuit.measure(my_alg.q_reg[1], my_alg.c_reg[1]) # measures second qubit\n",
    "\n",
    "# print list of gates in the circuit\n",
    "print('List of gates:')\n",
    "for circuit in my_alg.q_circuit:\n",
    "    print(circuit.name)\n",
    "\n",
    "#Execute the quantum algorithm\n",
    "result = my_alg.Q_program.execute(my_alg.circ_name, backend=my_alg.backend, shots= my_alg.shots)\n",
    "\n",
    "#Show the results obtained from the quantum algorithm \n",
    "counts = result.get_counts(my_alg.circ_name)\n",
    "\n",
    "print('\\nThe measured outcomes of the circuits are:',counts)\n",
    "\n",
    "# credits to:  https://github.com/QISKit/qiskit-tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.4 Description of the algorithm\n",
    "\n",
    "\n",
    "Shor's algorithm exploits both classical and quantum computation. Classical operations are carried out in steps where an efficient classical algorithm exists, while quantum operations are used to find the periodicity of the function needed to factor $N_0$. Here we present an outline of the algorithm.\n",
    "\n",
    "<img src=\"figures/8/shor1.jpeg\"  width=\"400\">\n",
    "$$\\text{3. Quantum circuit for Shor's algorithm.}$$\n",
    "\n",
    "\n",
    "To run the algorithm, we need two quantum registers. One contains the order or period, called period register, and the other contains the results of the computation, called computational register. The size of both registers depends on the number $N_0$ to be factored. In particular, the period register should contain a number of qubits $n_p$ in the interval $2 \\text{log} N_0 \\leq n_p < 2 \\text{log} (2N_0) $ and the computational register should be large enough to be able to represent the number $N_0-1$, thus $n_q = \\text{log} (N_0-1)$.\n",
    "\n",
    "First, one needs to check if $N_0$ is even. If $N_0$ is even, one of the factor is $2$ and the other is $N_0/2$. \n",
    "\n",
    "\n",
    "\n",
    "If $N_0$ is odd, a base $a$ is picked randomly among the numbers from $0$ to $N_0-1$. Then, check if $a$ is a factor of $N_0$, by checking if gcd$\\left(a,N_0 \\right) \\neq 1$. If $a$ is has a common divisor with $N_0$, then one factor is given by gcd$\\left(a,N_0 \\right)$ and the other by $N_0/$gcd$\\left(a,N_0 \\right)$. If $a$ is co-prime with $N_0$ ( gcd$\\left(a,N_0 \\right) = 1$ ) then one needs to compute the function $a^x$ mod $N_0$, the modular exponentiation function (MEF), for $x=1,2,3,..Q-1$ where $N_0^2 \\leq Q < 2N_0^2$ and $Q=2^{n_p}$. \n",
    "\n",
    "At the beginning of the quantum algorithm used to compute the MEF we need two registers initialized to zero $\\lvert 00...0 \\rangle \\lvert 00...0 \\rangle $. The first register also known as the period register will stores all the possible values of the exponent $x$ by creating a uniform superposition of all possible bit strings through Hadamard gates on all qubits $\\frac{1}{\\sqrt{Q}} \\sum_{x=0}^{Q-1} \\lvert x \\rangle $, and the second register, called the computational register will store the results of the MEF $a^x$ mod $N_0$, $\\lvert a^x $ mod $ N_0 \\rangle$. Thus, after the first step, one has\n",
    "\n",
    "\\begin{equation}\n",
    "\t\\frac{1}{\\sqrt{Q}} \\sum_{x=0}^{Q-1} \\lvert x \\rangle \\lvert a^x mod N_0 \\rangle\n",
    "\\end{equation}\n",
    "\n",
    "Then, we apply the quantum Fourier transform (QFT) to the first register, so that $\\lvert x \\rangle \\rightarrow \\frac{1}{\\sqrt{Q}} \\sum_{s=0}^{Q-1} e^{\\frac{2\\pi i s x}{Q}}\\lvert s \\rangle$. As a result of the QFT, interference between all the possible states occurs and only the periodic ones survive. That is, if one measures the first register, one will see a value of $s$ such that $\\frac{sx}{Q}$ is an integer $d$ when $x$ a multiple of the period $r$. Which means that $\\frac{s}{Q}=\\frac{d}{r}$. By knowing the fraction $\\frac{s}{Q}$ one can find the value of $r$ through the continued fraction algorithm.\n",
    "Now, if $r$ is odd or $r=0$, the algorithm fails and one needs to restart by picking a different base $a$. If $r$ is even, one can factorize $a^r - 1$  mod $N_0$ into $\\left( a^{\\frac{r}{2}} - 1\\right)\\left( a^{\\frac{r}{2}} + 1\\right)$ mod $N_0$. The final step is to check if $a^{\\frac{r}{2}} + 1 \\, \\text{mod} \\, N_0 \\neq -1$. If that's true, then $\\text{gcd}\\left(a^{\\frac{r}{2}} + 1,N_0 \\right)$ will be the first factor and $\\text{gcd}\\left(a^{\\frac{r}{2}} - 1,N_0 \\right)$ the other.\n",
    "\n",
    "\n",
    "The execution of this version of the algorithm requires $n=log_2 \\left( N_0 \\right)$ qubits in the computational register to perform the modular exponentiation and another $2n$ qubits in the period register to perform the QFT. Thus, the algorithm requires a total number of $3n$ qubits."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.5 Example of Shor's factoring algorithm\n",
    "\n",
    "Let us see here an example of Shor's factoring algorithm for $N_0=21$. Since $21$ needs five bits to be represented in binary, we need at least $\\text{log}_2Q = 2*5$ qubits for factoring. \n",
    "\n",
    "<ol>\n",
    "    \n",
    "    <li> Check if $N_0=21$ is even: $21 \\, \\text{mod} \\, 2=1$, $N_0$ is not even.</li>\n",
    "    \n",
    "    <li> Pick a base $a$ at random. Let's say $a=2$. </li>\n",
    "    \n",
    "    <li> Check if $a$ has any common factors with $N_0$: $\\text{gcd}\\left(a,N_0 \\right) = 1$, it doesn't.</li>\n",
    "    \n",
    "    <li> Initialize two qubit registers with $n_p = 2\\times5=10$ qubits, as $5$ qubits are needed to represent $21$ in binary. We call the first register the $period$ register and the second register the $computational$ register\n",
    "    $$ \\lvert \\psi \\rangle = \\lvert 0 \\rangle^{\\otimes 10}_p \\lvert 0 \\rangle^{\\otimes 10}_c $$ </li>\n",
    "    \n",
    "    <li> Apply Hadamard gate on all ther qubits of the first register to create a uniform superposition of all $2^{10}$ possible values \n",
    "    $$ \\lvert \\psi \\rangle =\\frac{1}{\\sqrt{2^{10}}} \\sum_{x=0}^{2^{10}-1} \\lvert x \\rangle_p \\lvert 0 \\rangle_c $$ </li>\n",
    "    \n",
    "    <li> Apply the modular exponentiation function $a^x \\, \\text{mod} \\, N_0$ on the second register, for each of the stored value of $x$ in the first register.\n",
    "        \n",
    "$$ \\lvert \\psi \\rangle =\\frac{1}{\\sqrt{1024}} \\sum_{x=0}^{1023} \\lvert x \\rangle_p \\lvert 2^x \\, \\text{mod} \\, 21 \\rangle_c\n",
    "        = \\\\ = \\frac{1}{\\sqrt{1024}} \\left( \\lvert 0 \\rangle_p \\lvert 1 \\rangle_c  + \\lvert 1 \\rangle_p \\lvert 2 \\rangle_c + \\lvert 2 \\rangle_p \\lvert 4 \\rangle_c + \\lvert 3 \\rangle_p \\lvert 8 \\rangle_c + \\lvert 4 \\rangle_p \\lvert 16 \\rangle_c  + \\lvert 5 \\rangle_p \\lvert 11 \\rangle_c + \\lvert 6 \\rangle_p \\lvert 1 \\rangle_c + \\lvert 7 \\rangle_p \\lvert 2 \\rangle_c + \\lvert 8 \\rangle_p \\lvert 4 \\rangle_c ... \\right) $$\n",
    "        By looking at the values stored in the second register, we can find out what is their periodicity. In particular, it can be seen that the values start repeating with order $r=6$ (value of the first register). However, we need to do a few more steps to allow the quantum computer to find the answer by itself.</li>\n",
    "     \n",
    "     <li> To simplify the example, we will adopt the \"Principle of implicit measurement\": Without loss of generality, any qubits which are not measured at the end of the quantum circuit may be assumed to be measured. Thus, let us use the principle of implicit measurement on the second register. Since each term of the superposition has equal weight, each outcome is equally likely, therefore one will see one of the following values: $ \\lvert 1 \\rangle_c$ , $\\lvert 2 \\rangle_c$, $ \\lvert 4 \\rangle_c$, $\\lvert 8 \\rangle_c$, $\\lvert 16 \\rangle_c$, $\\lvert 11 \\rangle_c$; with probability $1/6$.\n",
    "     Assume that the state $ \\lvert 4 \\rangle_c$ is measured, then we are left with the composite state:\n",
    "     \n",
    "     $$ \\lvert \\psi \\rangle = \\frac{\\sqrt{6}}{\\sqrt{1024}} \\left( \\lvert 2 \\rangle_p \\lvert 4 \\rangle_c + \\lvert 8 \\rangle_p \\lvert 4 \\rangle_c + \\lvert 14 \\rangle_p \\lvert 4 \\rangle_c + \\lvert 20 \\rangle_p \\lvert 4 \\rangle_c  ... \\right) $$</li>\n",
    "     \n",
    "     <li> Apply the QFT to the first register\n",
    "     $$ \\lvert x \\rangle_p \\rightarrow \\frac{1}{\\sqrt{1024}} \\sum_{s=0}^{1023} e^{2\\pi i \\frac{s x}{1024}}\\lvert s \\rangle_p $$\n",
    "     Which transform each of the terms in the period register as:\n",
    "     $$ \\lvert 2 \\rangle_p \\rightarrow \\frac{1}{\\sqrt{1024}} \\sum_{s=0}^{1023} e^{2\\pi i \\frac{2 s}{1024}}\\lvert s \\rangle_p $$\n",
    "     $$ \\lvert 8 \\rangle_p \\rightarrow \\frac{1}{\\sqrt{1024}} \\sum_{s=0}^{1023} e^{2\\pi i \\frac{8 s}{1024}}\\lvert s \\rangle_p $$\n",
    "     $$ \\lvert 14 \\rangle_p \\rightarrow \\frac{1}{\\sqrt{1024}} \\sum_{s=0}^{1023} e^{2\\pi i \\frac{14 s}{1024}}\\lvert s \\rangle_p $$\n",
    "    $$ \\lvert 20 \\rangle_p \\rightarrow \\frac{1}{\\sqrt{1024}} \\sum_{s=0}^{1023} e^{2\\pi i \\frac{20 s}{1024}}\\lvert s \\rangle_p $$\n",
    "    $$...$$\n",
    "    \n",
    "    So, we can write:\n",
    "    $$ \\lvert \\psi \\rangle = \\frac{\\sqrt{6}}{1024} \\sum_{s=0}^{1023} \\left( e^{2\\pi i \\frac{2s}{1024}} + e^{2\\pi i \\frac{8s}{1024}} + e^{2\\pi i \\frac{14s}{1024}} + e^{2\\pi i \\frac{20s}{1024}}  ... \\right)  \\lvert s \\rangle_p \\lvert 4 \\rangle_c  = \\\\\n",
    "    = \\frac{\\sqrt{6}}{1024} \\sum_{s=0}^{1023} e^{2\\pi i \\frac{2s}{1024}} \\left( 1 + e^{2\\pi i \\frac{6s}{1024}} + e^{2\\pi i \\frac{12s}{1024}} + e^{2\\pi i \\frac{18s}{1024}}  ... \\right)  \\lvert s \\rangle_p \\lvert 4 \\rangle_c$$\n",
    "    </li>\n",
    "\n",
    "    <li> Measure the period register. The probability of finding a certain value $\\lvert s \\rangle_p$ is:\n",
    "    $$ P(s) = \\lvert \\frac{\\sqrt{6}}{1024} e^{2\\pi i \\frac{2s}{1024}} \\left( 1 + e^{2\\pi i \\frac{6s}{1024}} + e^{2\\pi i \\frac{12s}{1024}} + e^{2\\pi i \\frac{18s}{1024}}  ... \\right) \\rvert^2 $$\n",
    "    \n",
    "    Because of the possible sign difference between all different terms, the values of $s$ which have the highest likelyhood to be observed are the ones for which the phase terms all have the same sign and they add up. That is, $\\frac{s}{1024}=\\frac{d}{6}$, where $d$ is an integer. Therefore, a value $s=\\frac{1024 \\cdot d}{6}$ where $d=1,2,3,...$ will most likely be observed.\n",
    "    Let us assume that the value $s=853$ is measured.\n",
    "    </li>\n",
    "\n",
    "    <li> To find the period from the value of $s$ measured, one then uses the continued fraction algorithm in the following way: we find the fraction $\\frac{d}{r}$ which approximates the fraction $\\frac{s}{Q}$ to a fixed precision $2*Q$\n",
    "    $$\\lvert \\frac{853}{1024} - \\frac{d}{r} \\rvert < \\frac{1}{2048} $$\n",
    "    So, let's find $d$ and $r$ with the continued fraction algorithm:\n",
    "    $$ \\frac{853}{1024} = 0 + \\frac{1}{1 + \\frac{1}{4+ \\frac{1}{1+\\frac{1}{1+\\frac{1}{1+\\frac{1}{84+\\frac{1}{2}}}}}}} $$\n",
    "   which gives as possible fractions $\\frac{d}{r}$: $\\frac{1}{1}$, $\\frac{4}{5}$, $\\frac{5}{6}$, ...\n",
    "   \n",
    "   imposing the condition written above, we find that the only fraction satisfying it is $\\frac{d}{r} = \\frac{5}{6}$. Therefore, the period is $r=6$!\n",
    "    </li>\n",
    "    \n",
    "    <li> Once the period has been found, we can find the factors of $N_0=21$ almost immediately. First check if $r$ is even then check if $a^{\\frac{r}{2}} + 1 \\, \\text{mod} \\, N_0 \\neq -1 $ :\n",
    "    $$ 6 \\, \\text{mod} \\, 2 = 0 $$\n",
    "    $$  2^{\\frac{6}{2}} + 1 \\, \\text{mod} \\, 21 = 9$$\n",
    "    Therefore the two factors are: $ p = \\text{gcd}\\left(a^{\\frac{r}{2}} + 1,N_0 \\right)$ and $q = \\text{gcd}\\left(a^{\\frac{r}{2}} - 1,N_0 \\right)$\n",
    "    $$p = \\text{gcd}\\left(9,21 \\right) = 3$$\n",
    "    $$q = \\text{gcd}\\left(7,21 \\right) = 7$$\n",
    "    </li>\n",
    "\n",
    "</ol>\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "\n",
    "<ol>\n",
    "<li>\n",
    "Calculate the following modular expression:\n",
    "\n",
    "<ol>\n",
    "<li>\n",
    "$ 54 \\, \\text{(mod 5)}$\n",
    "</li>\n",
    "<li>\n",
    "$ 54 \\times 35 \\, \\text{(mod 17)}$\n",
    "</li>\n",
    "<li>\n",
    "$ 9^{2} \\, \\text{(mod 6 )}$\n",
    "</li>\n",
    "</ol>\n",
    "</li>\n",
    "\n",
    "<li>\n",
    "Rewrite the following fractions using the continued fraction algorithm\n",
    "\n",
    "<ol>\n",
    "<li>\n",
    "$ \\frac{45}{251}$\n",
    "</li>\n",
    "<li>\n",
    "$ \\frac{71}{385}$\n",
    "</li>\n",
    "<li>\n",
    "$ \\frac{512}{1027}$\n",
    "</li>\n",
    "</ol>\n",
    "\n",
    "</li>\n",
    "\n",
    "<li>\n",
    "What is the minimum number of qubits needed to factor:\n",
    "\n",
    "<ol>\n",
    "<li>\n",
    "$ 15 $\n",
    "</li>\n",
    "\n",
    "<li>\n",
    "$ 4,367,398 $\n",
    "</li>\n",
    "</ol>\n",
    "\n",
    "</li>\n",
    "\n",
    "<li>\n",
    "Derive the expression of the Quantum Fourier Transform of three qubits, similarly to how it is done in the example shown in Section 8.3\n",
    "\n",
    "</li>\n",
    "\n",
    "<li>\n",
    "Following the steps highlighted in the example given in Section 8.5, show the factoring of $N_0=15$ with $a=2$.\n",
    "</li>\n",
    "\n",
    "<li>\n",
    "Write a QISKit program that calculates the QFT of three qubits in the state $\\lvert 010 \\rangle$ .\n",
    "</li>\n",
    "\n",
    "\n",
    "</ol>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "\n",
    "[1] P. Shor, in Proc. 35th Annu. Symp. on the Foundations of\n",
    "Computer Science, edited by S. Goldwasser (IEEE Computer\n",
    "Society Press, Los Alamitos, California, 1994),\n",
    "p. 124-134.\n",
    "\n",
    "\n",
    "[2] L. M. K. Vandersypen, M. Steffen, G. Breyta, C. S. Yannoni,\n",
    "M. H. Sherwood, and I. L. Chuang, Nature 414,\n",
    "883 (2001).\n",
    "\n",
    "[3] Monz, T. et al., Science 351, 1068-1070 (2016).\n",
    "\n",
    "\n",
    "[4] C. Y. Lu, D. E. Browne, T. Yang, and J. W. Pan, Physical\n",
    "Review Letters 99, 250504 (2007).\n",
    "\n",
    "[5] E. Martin-Lopez, A. Laing, T. Lawson, R. Alvarez, X.-Q.\n",
    "Zhou, and J. L. O'Brien, Nat Photon advance online publication (2012).\n",
    "\n",
    "[6] Lanyon, B. P. et al., Phys. Rev. Lett. 99, 250505 (2007).\n",
    "\n",
    "[7] A. Politi, J. C. F. Matthews, and J. L. O' Brien, Science 25, 1221 (2009).\n",
    "\n",
    "[8] E. Lucero, R. Barends, Y. Chen, J. Kelly, M. Mariantoni, A. Megrant, P. O'Malley, D. Sank, A. Vainsencher, J. Wenner, T. White, Y. Yin, A. N. Cleland and John M. Martinis, Nature Physics volume 8, 719-723 (2012).\n",
    "\n",
    "[9] P. J. Coles et al. (2018), arXiv:1804.03719.\n",
    "\n",
    "[10] J. A. Smolin, G. Smith, and A. Vargo, Nature (London) 499, 163-165 (2013).\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
